Collective non-Abelian instabilities in a melting Color Glass Condensate 
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We present first results for 3+1-D simulations of SU(2) Yang-Mills equations for matter expanding 
into the vacuum after a heavy ion collision. Violations of boost invariance cause a non-Abelian 
Weibel instability leading soft modes to grow with proper time r as exp(F^/ g 2 fir), where g 2 fu, is a 
scale arising from the saturation of gluons in the nuclear wavefunction. The scale for the growth 

rate F is set by a plasmon mass, defined as u p i — Ko \J generated dynamically in the collision. 

We compare the numerical ratio F/«o to the corresponding value predicted by the Hard Thermal 
Foop formalism for anisotropic plasmas. 



Experiments at the Relativistic Heavy Ion Collider 
(RHIC) at Brookhaven National Laboratory indicate 
that a thermalized Quark Gluon Plasma (QGP) may 
be formed in ultrarelativistic collisions of beams of gold 
ions . Phenomenological analyses of RHIC data point 
to thermalization of matter on early time scales of or- 
der 1 fm/c after the collision fj. Understanding how 
wavefunctions of high energy nuclei decohere completely 
on short time scales to form a thermalized QGP is an 
outstanding theoretical puzzle of great interest. 

At central rapidities, where thermalization is most 
likely, the wavefunctions of colliding nuclei are domi- 
nated by small x modes which, due to their large oc- 
cupation numbers, are described as classical fields 0. 
These fields, and their evolution with energy, can be com- 
puted in a classical effective field theory called the Color 
Glass Condensate (CGC) characterized by a semi- 
hard scale, Q s (x) » Aqcd- This "saturation" scale en- 
sures that the dynamics can in principle be understood 
in weak coupling; it is estimated [4| to be Q s w 1.4 GeV 
for RHIC energies and Q s w 2.2 GeV at Large Hadron 
Collider (LHC) energies of y/s — 5.5 TeV/nucleon. The 
initial conditions for the collision can be expressed self- 
consistently in terms of the computed classical fields of 
each of the nuclei [5j . 

The classical dynamics of fields produced in the col- 
lision is approximately boost invariant. Assuming it to 
be strictly so, the resulting 2+1 dimensional Yang-Mills 
equations generating the space-time dynamics of fields 
(the "melting of the CGC" ) was investigated numerically, 
and the energy and number distributions of the classically 
produced gluons computed On proper time scales 
r ~ 1/Qs, the energy density behaves as s ~ 1/r, sug- 
gesting free streaming of gluons in the transverse plane at 
these early times. In their "bottom up" scenario, Baicr 
et al. suggested that thermalization is a consequence of 
subsequent re-scattering of on-shell gluons by 2 <-> 2 Q 
and 2 <-* 3 processes ja The bottom up estimates 
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This thermalization time scale, 



at RHIC energies, is significantly larger than 1 fm/c. 

Scattering rates in the bottom up scenario are con- 
trolled by the Debye mass, which sets the range of par- 
ticle interactions. Interestingly, the Debye mass squared 
can change sign for anisotropic "CGC like" initial con- 
ditions, giving rise to an imaginary component in the 
gluon dispersion relation |14| . This corresponds to the 
non-Abelian analog of the collective Weibel instabil- 
ity [loL HH Il2^ leading to an exponential growth of soft 
modes. This phenomenon may invalidate or at least mod- 
ify [13j bottom up thermalization. Recently, several nu- 
merical studies have explored the behavior of instabili- 
ties in the finite temperature Hard Thermal Loop (HTL) 
formalism extended to anisotropic non-Abelian plas- 
mas These studies find exponentially growing soft 
modes. However, the most recent detailed 3+1-D simu- 
lations find that the exponential modes saturate at late 
times when non-linear self-interactions become impor- 
tant. Also interesting is a subsequent phenomenon, anal- 
ogous to Kolmogorov cascading in turbulent plasmas, 
where energy is transferred back to hard modes |l7| . 
While all these studies consider anisotropic particle dis- 
tributions a la the CGC, neither the expansion of the 
system nor the specific CGC dynamics following from 
physics of the small x modes is implemented. 

We present in this letter, for an SU(2) Yang-Mills 
theory, a first 3+1-dimensional numerical study of non- 
Abelian collective instabilities generated in the expand- 
ing CGC. The three spatial dimensions here are the 
two transverse co-ordinates x±, the space time rapidity 
i] = |m(|^). The fourth dimension is given by the 

proper time r = s/t 2 — z 2 . The strict boost invariance 
of fields (their lack of dependence on 77) in Ref. [|| is 
relaxed [2Cf. Violations of boost invariance arise from 
finite energy constraints and from small x quantum cor- 
rections. The magnitude of the violations from the for- 
mer are small. The latter are significantly larger and are 



of order unity over rapidity intervals of order Y ~ 
We will only consider small violations of boost invari- 
ance here; as a consequence, the gluon field configura- 
tions arc highly anisotro pic in momentum space. Ki- 
netic theory studies [lH Il2l Hif show that anisotropic 
distributions of hard modes rather than details of their 
dynamics drive the non-Abelian Weibel instability. Our 
initial conditions are therefore similar, in this respect, 
because the high momentum anisotropically distributed 
hard modes of the field with k± ~ Q s play the role of 
"hard" particles in their coupling to the soft modes with 
k± w T~ l k v << Q s . The effects of large amplitude vio- 
lations of boost invariance will be commented on briefly 
and discussed further in a follow up to this work 18]. 

In A T = gauge, initial conditions that violate strict 
boost invariance for the dynamical fields A± , A v and 
their canonically conjugate momenta Ex, E v are 

Ai = A ; A v = ; E. L = SE. L ; E n = £ n + SE n , 

where, for each configuration of the color charges of the 
two nuclei, Ai = Ai(x±), A n = 0, Ei — 0, £ v = £ v (x±) 
are the boost invariant initial conditions studied in pre- 
vious simulations 0. The rapidity dependent functions 
SEi and SE V are constructed to satisfy Gauss's law, 
DiSEi + D V E V = 0, at the initial time r = tq. For each 
SEi, random configurations 6Ei(x±) are drawn, satisfy- 
ing (SEi(x±_) SEi(y±)) = S^(x±_ — y±). These random 
configurations are multiplied by another random function 
f(rj) = d v F(rj). One obtains, at t = tq, 

5Ei(x ± ,r)) = f{ri) 5Ei(x±) ; 5E V = -F(rj) A SE^xx) . 

The fields F(rj) are Gaussian white noise distributed, 
(F{rj) F{r}')) — A 2 S(r] — ?y'), with the amplitude of vi- 
olations of boost invariance governed by the parameter 
A. These initial conditions are not unique. Their virtue 
is that Gauss' law is manifest and periodic boundary con- 
ditions can be applied in the 77 direction. A broader class 
of initial conditions will be considered in future. 

Our numerical procedure is as follows |2l|. The small 
x classical fields before the collision are determined from 
their respective classical color source densities by solv- 
ing Poisson's equations Q. These color charge distribu- 
tions are Gaussian distributed with a variance g 2 [i. This 
momentum scale (closely related and similar in magni- 
tude to Q s ) and the size of the system L are the 
physical dimensionful scales in the problem. For periodic 
boundary conditions, L 2 = irR 2 , where R is the radius 
of the nuclei. For each configuration of color charges, the 
initial conditions are determined, and an adaptive leap- 
frog algorithm evolves Hamilton's equations for the dy- 
namical fields and their canonically conjugate momenta 
on a discretized space-time lattice. Physical results, ob- 
tained by averaging the results over all color configura- 
tions of the sources, are expressed in terms of g 2 [i and the 
dimensionless combination g 2 fiL. Higher energies and/or 
larger nuclei correspond to larger values of g 2 [iL. 
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FIG. 1: The maximum Fourier mode amplitudes of t 2 T vv 
for g 2 fiL = 67.9, N± = N n = 64, N v a v = 1.6. Also shown 
are best fits with exp r and exp ^fr behavior. The former is 
clearly ruled out by the data. 

The lattice parameters, in dimensionless units, are a) 
N±_ and , the number of lattice sites in the x± and r) 
directions respectively; b) g 2 fia± and a n , the respective 
lattice spacings; c) To/a± and St, the time at which the 
simulations are initiated and the stepping size respec- 
tively; and finally d) A, the initial size of the rapidity 
fluctuations. The continuum limit is obtained by hold- 
ing the physical combinations g 2 \x a± N± = g 2 /iL and 
a v N n — L v fixed while sending St, g 2 fi a± and a n to 
zero. For this study, we pick L v = 1.6 units of rapidity. 
Variations of L n will be commented on later. The magni- 
tude of violations of boost invariance, as represented by 
A, is a physical quantity; here, we study results for small 
values of A ~ 1CP 11 — 10~ 8 . The initial time is chosen to 
ensure that for A = 0, we recover earlier 2+1-D results; 
we set To = 0.05 a±. Our results are insensitive to vari- 
ations that are a factor of 2 larger or smaller than this 
choice. 

To study the growth rate of instabilities due to viola- 
tions of boost invariance, we define 

f fci - 0) - J dr, e X p(ir ) k ri ){T^(xx,v))±, P , 

where denote components of the stress energy ten- 
sor and ()_L, P denotes an average over the transverse co- 
ordinates and over all color charge configurations respec- 
tively. We will look in particular at this Fourier transform 
for the longitudinal pressure r 2 T nn . The magnitude of 
this quantity is a useful measure of isotropization; mo- 
mentum distributions are isotropic when 2 T 2 T riri / [T xx + 
T vv ) -> 1. In Fig.Q] we plot the maximal value of t 2 X"™ 
at each time step, as a function of g 2 ! /jt. The data are 
for a 64 3 lattice and correspond to g 2 /J,L — 67.9 and 
L, ; = 1.6. The maximal value remains nearly constant 
until g 2 [it w 250, beyond which it grows rapidly. A 
best fit to the functional form cq + ciexp(c2T C3 ) gives 
C2 = 0.427 ± 0.01 for C3 = 0.5; the coefficients cq, c\ are 




FIG. 2: The mode number k v for which the maximum am- FIG. 3: Time evolution of u> pU for fixed g 2 /xL = 22.6 and lat- 

plitude of T 2 T vn occurs (see Fig[TJ; again for g 2 (j,L = 67.9, tice spacings g 2 fia ± = 0.707,0.354,0.177 (N± = 32,64, 128), 

N± = N v — 64, N v a v — 1.6. fc m i n denotes the smallest k n respectively, 
mode that could be excited in the lattice simulation. 



small numbers proportional to the initial seed. It is clear 
from Fig. ^that the form exp(r \/g 2 ~Jvr) is preferred to a 
fit with an exponential growth in r. 

We digress briefly to note that the relevant time scales 
in Fig. n f° r RHIC collisions where <? 2 /i ~ 1 GeV, are 
too large to be interest. These however correspond to the 
time necessary for extremely small violations of boost 
invariance to become visible. The precise behavior of 
non-Abelian Weibel instabilities in a melting CGC are 
more transparent for small initial seeds. 

In Fig. |2 we plot the most unstable k v mode, corre- 
sponding to the maximal r 2 T m in Fig. ^ as a function 
of g 2 /J,r. Our result suggests that the soft k v modes are 
sensitive to the non-Abelian Weibel instability while the 
hard modes are not. Further, the maximal mode num- 
ber k v grows as a function of time. The amplitude of the 
soft modes dominates the spectrum for g 2 /J,r > 250. This 
timescale, however, is weakly dependent on our choice of 
L, ; = 1.6; raising L v to L v = 3.2 lowers the timescale 
for the onset of growth by some 20%, while the extracted 
growth rate stays the same. The large volume limit will 
be examined further in an upcoming work |l8|. 

We now turn to the correspondence between our re- 
sults and the HTL formalism for anisotropic plasmas. 
In the latter, the maximal Fourier amplitudes of compo- 
nents of the stress-energy tensor grow as exp(27r), where 



the growth rate 7 satisfies the relation 7 = J i 
maximally anisotropic particle distributions [12j. 
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where / is the anisotropic single particle distribution of 
the hard modes. For both isotropic as well as anistropic 
distributions (within the model used in jj^j, j2^|) m 2 ^ = 



5 W . 2 



pi- 



Interestingly, in close analogy to the HTL case, the 
exp(r^/ g 2 fJ.r) growth of the unstable soft modes in our 



case is closely related to a mass gap generated by the 
highly non-linear dynamics of soft modes in the expand- 
ing system. As in Ref. Q , fixing the residual gauge free- 
dom Vj_ • A±, the gluon dispersion is given by relation 




w(kj_) = - 



1 Tr[K i (k ±J B 4 (-k ± )+r 2 ^ t) (k ± )^ t) (-k ± )] 



r y Tr [^(kx)^(-kj.) + r- 2 A,(kj.)A,(-k±)] 

Remarkably, a mass-gap, which we associate with the 
plasmon mass lo p \ = uj(h± = 0), is generated. After 
initial transient behavior, it satisfies the relation 



^pi(r) = K g 2 fj, 



1 



The time evolution of this plasmon mass in units of 
w p i(r / 'g 2 ^i) 1 / 2 is shown in Fig. El for g 2 [iL = 22.6, with 
k = 0.3 ± 0.01; it is robust as one approaches the con- 
tinuum limit. The dependence of w p i on g 2 \i L is shown 
in Fig. HHil. 

We now assume the growth rate in the expanding sys- 
tem is modified as j s t a t T li T ) T i an d further assume, a) 
7( r ) = "ioo(V)/\/2, and b) to? x) (t) = fa^r), as in the 
static case. Since the plasmon mass is determined inde- 
pendently, these HTL relations predict the Fourier ampli- 
tude grows as exp (Vth. \J 9 ,2 M T J i with ^th. = V^^o- This 
HTL motivated value of the growth rate is compared to 
the growth rate extracted directly from our fits in the ta- 
ble for several values of g 2 [iL (computed for N± = 32, 64 
and 128, respectively). Remarkably, we find an agree- 
ment to within 4 — 6% accuracy. However, a consistent 
treatment gives the growth rate in the expanding case be 
exp(2 L dr'7(r')). If the HTL relation of 7 to the plas- 
mon mass is unchanged, an additional factor of 2 is ob- 
tained in T t h. relative to Tfn- It is not obvious that HTL 
relations generalize to the case of CGC initial conditions. 
In particular, studies are in progress to investigate how 
our measured uj p i relates to the HTL plasmon frequency. 
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FIG. 4: Dependence of uj p \, with respect to g 2 fJ,L; g 2 fia± — 



0.707 for all g 2 fiL except g 2 fJ.L 
1.06. 



67.9, for which it is g /xaj 
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Further refinements such as the approach to the con- 
tinuum limit in the 77 direction require larger lattices. 
They are especially important for studies where the 
scale A, governing violations of boost invariance, is 
increased. Much larger values of A are required to 
study whether isotropization of particle distributions oc- 
curs through classical collective instabilities, a combined 
instability/modified bottom up scenario |l3| or non- 
perturbative scenarios [19j. 

To summarize, we performed 3±1-D numerical sim- 
ulations of SU(2) Yang-Mills equations with CGC ini- 
tial conditions that describe the pre-equilibrium stage 
of high energy heavy ion collisions. Violations of strict 
boost invariance cause a non-Abelian Weibel instabil- 
ity. We demonstrated unambiguously that amplitudes 
of soft modes grow as exp(T^/ g 2 ^r) for gluonic matter 
expanding into the vacuum at the speed of light. This 
behavior differs from the exp(^ t) growth of instabilities 
in a stationary plasma. The scale governing the growth 
rate is determined by a plasmon mass generated by the 
non-Abelian dynamics. Extrapolating HTL motivated 
assumptions to an expanding system, we compared the 
growth rate determined from the plasmon mass to the 
growth rate extracted directly in our simulations. An- 
alytical and numerical work is in progress to test these 
assumptions. The possible isotropization of distributions 
due to collective instabilities arising from larger viola- 
tions of boost invariance is also under investigation. 

RV's research is supported by DOE Contract No. DE- 
AC02-98CH10886. He thanks the Alexander von Hum- 
boldt foundation for support during the early stages of 



this work. PR was supported by DFG-Forschergruppe 
EN164/4-4. We thank D. Bodeker, J. Engels, F. Gelis, 
D. Kharzeev, A. Krasnitz, M. Laine, T. Lappi, L. McLer- 
ran, Y. Nara and M. Strickland for fruitful discussions. 



[1] K. Adcox et al., [PHENIX collaboration], Nucl. Phys. 
A757 (2005) 184-283; J. Adams et al., [STAR collabora- 
tion], ibid., 102-183; B. B. Back et al., [PHOBOS collab- 
oration], ibid., 28-101; I. Arsene et al., [BRAHMS collab- 
oration], ibid ., 1-27. 

P. Huovinen, arXiv:nucl-th/0305064 in Hwa, R.C. (ed.) 
et al.: in QGP3 (World Scientific Publishers). 
L. McLerran, R. Venugopalan, Phys. Rev. D 49, 2233, 
3352 (1994); ibid 50, 2225 (1994). 

E. Iancu, R. Venugopalan, hep-ph/0303204 in Hwa, R.C. 
(ed.) et al.: QGP3 (World Scientific Publishers). 
A. Kovner, L. D. McLerran and H. Weigert, Phys. Rev. 
D 52, 3809 (1995); ibid., 6231. 

A. Krasnitz and R. Venugopalan, Nucl. Phys. B 557, 
237 (1999); Phys. Rev. Lett. 84, 4309 (2000); 86, 1717 
(2001); A. Krasnitz, Y. Nara and R. Venugopalan, Phys. 
Rev. Lett. 87, 192302 (2001); Nucl. Phys. A 717, 268 
(2003); 727, 427 (2003); T L appi, Phys. Rev. C 67, 
054903 (2003); arXiv:hep-ph/0505095 
A. H. Mueller, Nucl. Phys. B 572, 227 (2000); Phys. 
Lett. B 475, 220 (2000); J. Bjoraker and R. Venugopalan, 
Phys. Rev. C 63, 024609 (2001). 

S. M. H. Wong, Phys. Rev. C 54, 2588 (1996); Z. Xu and 
C. Greiner, Phys. Rev. C 71, 064901 (2005). 
R. Baier, A.H. Mueller, D. Schiff, D.T. Son, Phys. Lett. 
B 502, 51 (2001). 

E.S. Weibel, Phys. Rev. Lett 2, 83 (1959) 
S. Mrowczynski, Phys. Lett. B 214, 587 (1988); Phys. 
Lett. B 314, 118 (1993); Phys. Lett. B 363, 26 (1997); 
J. Randrup, S. Mrowczynski, Phys. Rev. C 68, 034909 
(2003). 

P. Arnold, J. Lenaghan, G.D. Moore, JHEP 0308, 
002 (2003); P. Arnold, J. Lenaghan, G. D. Moore and 
L. G. Yaffe, Phys. Rev. Lett. 94, 072302 (2005). 
A.H. Mueller, A.I. Shoshi S.M.H. Wong, arX iv: 
hep-ph/0505164 D. Bodeker, arXiv:hep-ph/0508223 
P. Romatschke, M. Strickland, Phys. Rev. D 68, 036004 
(2003); Phys. Rev. D 70, 116006 (2004). 
E. Braaten, R.D. Pisarski, Phys.Rev.Lett. 64, 1338 
(1990). 

A. Rebhan, P. Romatschke, M. Strickland, Phys. Rev. 
Lett 94, 102303 (2005); JHEP 0509, 041 (2005). P. 
Arnold, G.D. Moore, L.G. Yaffe, Phys.Rev.D 72, 054003 
(2005); for a related approach, see A. Dumitru and 

Y. Nara, Phys. Lett. B 621, 89 ( 2005). 

P. Arnold and G. D. Moore, arXiv:hep-ph/0509206 
hep-ph/0509226 

P. Romatschke and R. Venugopalan, in preparation. 
Y. V. K ovchegov, arXiv:hep-ph/0507134 
hep-ph/0503038|| D. Kharzeev and K. Tuchin, Nucl' 
Phys. A 753. 316 (2005). 

It is essential to distinguish "strict" boost invariance of 
fields from that of single particle distributions. The latter 
can be preserved even if the former is violated. 
This procedure is similar to that of Ref. |(|. Further de- 
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tails will be presented in Ref. [l8| . 
[22] In the classical theory of the CGC, Q 2 = g 4 fi 2 N c /2iv ■ 

ln(g 2 fi/A), where A is an infrared scale of order Aqcd- 
[23] Averaging over longitudinal and transverse modes (see 

the 2nd paper by Rebhan et al. in Ref. [Ig), one obtains 

w pl — ^jarctany^ and m 2 ^ = ^^arctany 7 ^, where £ 

denotes the strength of the anisotropy and m 2 is a soft 
scale proportional to the Debye mass. Therefore, lo^ = 
2m 2 



[24] The effect of small longitudinal fluctuations on transverse 
quantities should be rather small. We therefore extract 
u>(k±), for computational convenience, from a 2+1-D 
simulation. 

[25] The mass gap was noted in Ref. [g in context of infrared 
finite number distributions. Interestingly, the behavior of 
ui p i versus g 2 fJ,L is similar to that of the gluon liberation 
factor /jv. 



